Fisk distribution (fisk) — the log-logistic workhorse

The Fisk distribution (SciPy: scipy.stats.fisk) is widely known in applied statistics as the log-logistic distribution. It is a continuous distribution on the positive real line with Pareto-like tails.

A key representation makes much of its behavior intuitive:

If (X \sim \mathrm{Fisk}(c, \text{scale}=s)) with loc=0, then

[ \log X \sim \mathrm{Logistic}(\mu=\log s,; \sigma = 1/c). ]

So you can think of Fisk as “a logistic distribution on the log-scale”.

What you’ll learn

• how to write the PDF/CDF/quantile function (and how SciPy parameterizes them)

• which moments exist (and why many do not)

• how to sample with NumPy only via inverse-CDF / logistic tricks

• how to fit and validate the model with scipy.stats.fisk

---

import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import optimize, stats

pio.templates.default = "plotly_white"
# CKC convention for Plotly notebooks
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)

# Example parameters used throughout
c0 = 5.0
scale0 = 2.0
loc0 = 0.0

Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0

1) Title & Classification

• Name: fisk (Fisk / log-logistic distribution; SciPy: scipy.stats.fisk)

• Type: Continuous

• Support (standard): (x \in (0, \infty))

• Parameter space (standard): shape (c>0), scale (s>0)

• SciPy location/scale: loc \in \mathbb{R}, scale > 0 with [ X = \mathrm{loc} + \mathrm{scale},Z,\quad Z \sim \mathrm{Fisk}(c, 1). ]

We write (2-parameter form):

[ X \sim \mathrm{Fisk}(c, s). ]

The standard form is (\mathrm{Fisk}(c,1)), i.e. stats.fisk(c, loc=0, scale=1).

2) Intuition & Motivation

2.1 What it models

The Fisk/log-logistic distribution models strictly positive quantities where:

• values can vary over orders of magnitude, and

• the right tail is heavy (polynomial decay), and

• on the log-scale, the distribution looks approximately logistic.

A nice way to see the logistic link is to rewrite the CDF in terms of (y = \log x) (with loc=0):

[ F(x) = \frac{1}{1 + (s/x)^c} = \sigma\bigl(c(\log x - \log s)\bigr), \qquad \sigma(t) = \frac{1}{1+e^{-t}}. ]

So the “transition” happens around (x \approx s) (because (\log x - \log s \approx 0)).

2.2 Typical real-world use cases

• Survival analysis / reliability: time-to-event with a hazard that can rise and fall.

• Economics: income/wealth-like quantities with heavy tails.

• Hydrology / environmental extremes: positive heavy-tailed magnitudes.

• Generative modeling: as a simple positive heavy-tailed component.

2.3 Relations to other distributions

• Log-logistic: Fisk is the log-logistic distribution.

• Logistic on the log-scale: (\log X) is logistic.

• Burr Type XII: Fisk is a special case of Burr XII (shape parameter (k=1)).

• Pareto-like tails: (\mathbb{P}(X>x) \sim (s/x)^c) for large (x), so (c) behaves like a tail index.

3) Formal Definition

We use the common 2-parameter form with shape (c>0) and scale (s>0).

3.1 PDF

For (x>0):

[ f(x\mid c,s) = \frac{c}{s},\frac{(x/s)^{c-1}}{\bigl(1 + (x/s)^c\bigr)^2}. ]

3.2 CDF

For (x>0):

[ F(x\mid c,s) = \frac{1}{1 + (s/x)^c} = \frac{(x/s)^c}{1 + (x/s)^c} = \sigma\bigl(c(\log x - \log s)\bigr). ]

The survival function is

[ \bar F(x) = 1 - F(x) = \frac{1}{1 + (x/s)^c}. ]

3.3 Quantile function (inverse CDF)

For (u\in(0,1)):

[ Q(u) = F^{-1}(u) = s\left(\frac{u}{1-u}\right)^{1/c}. ]

3.4 SciPy parameterization

SciPy uses:

stats.fisk(c, loc=0, scale=1)

with (z = (x-\mathrm{loc})/\mathrm{scale}) and support (z>0).

def _log1pexp(a: np.ndarray) -> np.ndarray:
    """Stable log(1+exp(a)) for real a (NumPy-only)."""

    a = np.asarray(a, dtype=float)
    out = np.empty_like(a, dtype=float)

    pos = a > 0
    out[pos] = a[pos] + np.log1p(np.exp(-a[pos]))
    out[~pos] = np.log1p(np.exp(a[~pos]))
    return out


def _expit(a: np.ndarray) -> np.ndarray:
    """Stable logistic sigmoid 1/(1+exp(-a)) (NumPy-only)."""

    a = np.asarray(a, dtype=float)
    out = np.empty_like(a, dtype=float)

    pos = a >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-a[pos]))

    ea = np.exp(a[~pos])
    out[~pos] = ea / (1.0 + ea)
    return out


def _logit(p: np.ndarray) -> np.ndarray:
    """Stable log(p/(1-p)) for p in (0,1)."""

    p = np.asarray(p, dtype=float)
    return np.log(p) - np.log1p(-p)


def fisk_logpdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Log-PDF of Fisk(c, loc, scale) in SciPy's parameterization."""

    x = np.asarray(x, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale

    out = np.full_like(z, fill_value=-np.inf, dtype=float)
    mask = z > 0
    if not np.any(mask):
        return out

    logz = np.log(z[mask])
    a = c * logz  # log(z^c)

    out[mask] = (
        np.log(c)
        - np.log(scale)
        + (c - 1.0) * logz
        - 2.0 * _log1pexp(a)
    )

    return out


def fisk_pdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """PDF of Fisk(c, loc, scale)."""

    return np.exp(fisk_logpdf(x, c=c, loc=loc, scale=scale))


def fisk_cdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """CDF of Fisk(c, loc, scale)."""

    x = np.asarray(x, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale

    out = np.zeros_like(z, dtype=float)
    mask = z > 0
    if not np.any(mask):
        return out

    logz = np.log(z[mask])
    out[mask] = _expit(c * logz)
    return out


def fisk_sf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Survival function 1 - CDF of Fisk(c, loc, scale)."""

    x = np.asarray(x, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale

    out = np.ones_like(z, dtype=float)
    out[z <= 0] = 1.0

    mask = z > 0
    if not np.any(mask):
        return out

    logz = np.log(z[mask])
    out[mask] = _expit(-c * logz)
    return out


def fisk_ppf(u: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Quantile function (inverse CDF) of Fisk(c, loc, scale)."""

    u = np.asarray(u, dtype=float)

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    if np.any((u < 0) | (u > 1)):
        raise ValueError("u must be in [0, 1]")

    out = np.full_like(u, fill_value=np.nan, dtype=float)
    out[u == 0] = loc
    out[u == 1] = np.inf

    mask = (u > 0) & (u < 1)
    if not np.any(mask):
        return out

    out[mask] = loc + scale * np.exp(_logit(u[mask]) / c)
    return out

4) Moments & Properties

A useful fact is that the right tail is polynomial:

[ \mathbb{P}(X>x) = \bar F(x) = \frac{1}{1 + (x/s)^c} \sim (s/x)^c\quad (x\to\infty). ]

So (c) is a tail index: smaller (c) means heavier tails.

4.1 Raw moments (and existence)

For (k>0), the raw moment exists iff (c > k), and

[ \mathbb{E}[X^k] = s^k,\Gamma!\left(1+\frac{k}{c}\right)\Gamma!\left(1-\frac{k}{c}\right) = s^k,\frac{\frac{k\pi}{c}}{\sin\left(\frac{k\pi}{c}\right)}. ]

Consequences:

• Mean exists for (c>1): [\mathbb{E}[X] = s,\frac{\pi/c}{\sin(\pi/c)}.]

• Variance exists for (c>2): [\mathrm{Var}(X) = s^2\left(\frac{2\pi/c}{\sin(2\pi/c)} - \left(\frac{\pi/c}{\sin(\pi/c)}\right)^2\right).]

• Skewness exists for (c>3); kurtosis exists for (c>4).

4.2 Location summaries

• Median: (\mathrm{median}(X)=s) (and with SciPy loc, the median is loc + scale).

• Mode (for (c>1)): [\mathrm{mode}(X)=s\left(\frac{c-1}{c+1}\right)^{1/c}.] For (c\le 1), the density is decreasing and the mode is at the lower endpoint.

4.3 MGF / characteristic function

• The MGF (M_X(t)=\mathbb{E}[e^{tX}]) does not exist for any (t>0) (polynomial tails cannot beat exponential growth).

• The characteristic function (\varphi_X(t)=\mathbb{E}[e^{itX}]) exists, but it does not simplify to an elementary closed form; it’s typically computed numerically.

4.4 Entropy

Using the logistic-on-log-scale representation, the differential entropy is

[ h(X) = 2 + \log\left(\frac{s}{c}\right)\quad \text{(nats)}. ]

(Translation by loc does not change entropy; scaling multiplies it by adding (\log(\mathrm{scale})).)

def fisk_raw_moment(k: float, c: float, scale: float = 1.0) -> float:
    """Raw moment E[X^k] for Fisk(c, scale) (loc=0).

    Returns np.inf if the moment diverges (c <= k).
    """

    k = float(k)
    c = float(c)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")
    if k == 0:
        return 1.0

    if c <= k:
        return float(np.inf)

    a = np.pi * k / c
    return float((scale**k) * (a / np.sin(a)))


def fisk_entropy(c: float, scale: float = 1.0) -> float:
    """Differential entropy in nats for Fisk(c, scale) (loc drops out)."""

    c = float(c)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    return float(2.0 + np.log(scale / c))


def fisk_mode(c: float, loc: float = 0.0, scale: float = 1.0) -> float:
    """Mode of Fisk(c, loc, scale)."""

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    if c <= 1:
        return float(loc)  # density decreases from the boundary

    return float(loc + scale * ((c - 1.0) / (c + 1.0)) ** (1.0 / c))


def fisk_moments(c: float, loc: float = 0.0, scale: float = 1.0) -> dict:
    """Mean/variance/skewness/kurtosis (if they exist) + useful summaries."""

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    m1 = fisk_raw_moment(1.0, c=c, scale=scale)
    m2 = fisk_raw_moment(2.0, c=c, scale=scale)
    m3 = fisk_raw_moment(3.0, c=c, scale=scale)
    m4 = fisk_raw_moment(4.0, c=c, scale=scale)

    mean = loc + m1 if np.isfinite(m1) else np.inf

    if np.isfinite(m1) and np.isfinite(m2):
        var = m2 - m1**2
    else:
        var = np.inf

    skew = np.nan
    excess_kurt = np.nan

    if np.isfinite(m1) and np.isfinite(m2) and np.isfinite(m3) and var > 0:
        mu3 = m3 - 3 * m2 * m1 + 2 * (m1**3)
        skew = mu3 / (var ** 1.5)

    if np.isfinite(m1) and np.isfinite(m2) and np.isfinite(m3) and np.isfinite(m4) and var > 0:
        mu4 = m4 - 4 * m3 * m1 + 6 * m2 * (m1**2) - 3 * (m1**4)
        excess_kurt = mu4 / (var**2) - 3.0

    median = loc + scale

    return {
        "mean": float(mean),
        "var": float(var),
        "skew": float(skew) if np.isfinite(skew) else np.nan,
        "excess_kurtosis": float(excess_kurt) if np.isfinite(excess_kurt) else np.nan,
        "median": float(median),
        "mode": fisk_mode(c=c, loc=loc, scale=scale),
        "entropy": fisk_entropy(c=c, scale=scale),
    }


m0 = fisk_moments(c0, loc=loc0, scale=scale0)
m0

{'mean': 2.13791866423119,
 'var': 0.7145293838425228,
 'skew': 2.4852755496867136,
 'excess_kurtosis': 26.556191909249108,
 'median': 2.0,
 'mode': 1.8442158229634555,
 'entropy': 1.083709268125845}

5) Parameter Interpretation

5.1 Scale s

• scale = s sets the median: (\mathrm{median}(X)=s) (when loc=0).

• On the log-scale, it is a location parameter: (\mu = \log s).

5.2 Shape c

• c controls the steepness around the median and the tail heaviness.

• The tail index is c: (\mathbb{P}(X>x) \sim (s/x)^c).

• Moments exist only up to order < c.

• On the log-scale, c is the inverse of the logistic scale: (\sigma = 1/c).

5.3 Hazard shape (survival analysis intuition)

The hazard rate (h(x)=f(x)/(1-F(x))) (for loc=0) simplifies to

[ h(x) = \frac{c}{s},\frac{(x/s)^{c-1}}{1 + (x/s)^c}. ]

• For (c \le 1), the hazard is decreasing.

• For (c > 1), the hazard is typically unimodal (rises, then falls).

def fisk_hazard(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Hazard h(x) = f(x) / (1 - F(x))."""

    f = fisk_pdf(x, c=c, loc=loc, scale=scale)
    s = fisk_sf(x, c=c, loc=loc, scale=scale)

    out = np.full_like(f, fill_value=np.nan, dtype=float)
    mask = s > 0
    out[mask] = f[mask] / s[mask]
    return out


x = np.logspace(-3, 2, 800) * scale0
c_values = [0.7, 1.0, 2.0, 5.0]

fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=["PDF (log x)", "CDF (log x)", "Hazard (log x)"],
)

for c in c_values:
    fig.add_trace(
        go.Scatter(x=x, y=fisk_pdf(x, c=c, scale=scale0), name=f"c={c}", mode="lines"),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=x,
            y=fisk_cdf(x, c=c, scale=scale0),
            name=f"c={c}",
            mode="lines",
            showlegend=False,
        ),
        row=1,
        col=2,
    )
    fig.add_trace(
        go.Scatter(
            x=x,
            y=fisk_hazard(x, c=c, scale=scale0),
            name=f"c={c}",
            mode="lines",
            showlegend=False,
        ),
        row=1,
        col=3,
    )

for j in [1, 2, 3]:
    fig.update_xaxes(type="log", title_text="x", row=1, col=j)

fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_yaxes(title_text="h(x)", row=1, col=3)

fig.update_layout(title="Fisk distribution shape effects (scale fixed)", height=420)
fig.show()

6) Derivations

6.1 Expectation (raw moments)

Start from the definition for (k>0):

[ \mathbb{E}[X^k] = \int_0^\infty x^k,\frac{c}{s},\frac{(x/s)^{c-1}}{(1 + (x/s)^c)^2},dx. ]

Use the substitution (t = (x/s)^c), so (x = s,t^{1/c}) and (dx = \frac{s}{c} t^{1/c-1},dt). A pleasant simplification happens:

[ \frac{c}{s}(x/s)^{c-1},dx = \frac{c}{s},t^{(c-1)/c},\frac{s}{c}t^{1/c-1}dt = dt. ]

So

[ \mathbb{E}[X^k] = s^k \int_0^\infty \frac{t^{k/c}}{(1+t)^2},dt. ]

Recognize the Beta-function integral:

[ \int_0^\infty \frac{t^{p-1}}{(1+t)^{m}},dt = B(p, m-p),\quad 0<p<m. ]

Here (m=2) and (p=1+k/c), so the condition becomes (c>k). Then

[ \mathbb{E}[X^k] = s^k,B\left(1+\frac{k}{c},,1-\frac{k}{c}\right) = s^k,\Gamma\left(1+\frac{k}{c}\right)\Gamma\left(1-\frac{k}{c}\right). ]

Using (\Gamma(1+z)\Gamma(1-z)=\frac{\pi z}{\sin(\pi z)}) gives the sine form.

6.2 Variance

When (c>2),

[ \mathrm{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2, ]

with (\mathbb{E}[X]) requiring (c>1).

6.3 Likelihood (i.i.d. sample)

For i.i.d. observations (x_1,\dots,x_n > 0) (assuming loc=0), the log-likelihood is

[ \ell(c,s) = \sum_{i=1}^n \log f(x_i\mid c,s) = n\log c - n\log s + (c-1)\sum_i\log\left(\frac{x_i}{s}\right) - 2\sum_i \log\left(1+\left(\frac{x_i}{s}\right)^c\right). ]

There is no closed-form MLE for ((c,s)); we maximize (\ell) numerically. A common trick is to optimize over ((\log c, \log s)) to enforce positivity.

def fisk_loglikelihood(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> float:
    """Total log-likelihood for i.i.d. data under Fisk(c, loc, scale)."""

    return float(np.sum(fisk_logpdf(x, c=c, loc=loc, scale=scale)))


def fisk_mle_loc_fixed(x: np.ndarray, loc: float = 0.0) -> tuple[float, float]:
    """Numerical MLE for (c, scale) with loc fixed.

    Uses a logistic-on-log-scale initialization.
    """

    x = np.asarray(x, dtype=float)
    if np.any(x <= loc):
        raise ValueError("all observations must be > loc")

    x_pos = x - loc

    # Initialization from log-scale logistic approximation
    y = np.log(x_pos)
    mu_init = float(np.median(y))
    q25, q75 = np.quantile(y, [0.25, 0.75])
    iqr = float(q75 - q25)

    # Logistic IQR = 2*s*log(3), with s = 1/c
    c_init = float(max(2.0 * np.log(3.0) / max(iqr, 1e-6), 1e-3))
    scale_init = float(max(np.exp(mu_init), 1e-6))

    def nll(theta: np.ndarray) -> float:
        log_c, log_scale = float(theta[0]), float(theta[1])
        c = float(np.exp(log_c))
        scale = float(np.exp(log_scale))
        return -fisk_loglikelihood(x, c=c, loc=loc, scale=scale)

    res = optimize.minimize(
        nll,
        x0=np.array([np.log(c_init), np.log(scale_init)]),
        method="Nelder-Mead",
        options={"maxiter": 5000},
    )

    log_c_hat, log_scale_hat = res.x
    return float(np.exp(log_c_hat)), float(np.exp(log_scale_hat))


# Quick demo on synthetic data
x_demo = stats.fisk.rvs(c0, loc=loc0, scale=scale0, size=2000, random_state=rng)
c_hat, scale_hat = fisk_mle_loc_fixed(x_demo, loc=loc0)

print("true (c, scale) =", (c0, scale0))
print("MLE  (c, scale) =", (c_hat, scale_hat))

true (c, scale) = (5.0, 2.0)
MLE  (c, scale) = (5.081944016481647, 1.9992827793076853)

7) Sampling & Simulation

Inverse-transform sampling uses the quantile function.

If (U\sim\mathrm{Uniform}(0,1)), then

[ X = Q(U) = \mathrm{loc} + \mathrm{scale},\exp\left(\frac{\mathrm{logit}(U)}{c}\right) = \mathrm{loc} + \mathrm{scale}\left(\frac{U}{1-U}\right)^{1/c}. ]

Algorithm (NumPy-only):

1. Draw (u_1,\dots,u_n) i.i.d. from Uniform(0,1).

2. Clip away from exactly 0 or 1 (to avoid ±inf after logit).

3. Compute z = exp(logit(u) / c) and return loc + scale * z.

Because the distribution is logistic on the log-scale, this is equivalent to sampling

[ Y \sim \mathrm{Logistic}(\mu=\log(\mathrm{scale}), \sigma = 1/c),\qquad X = \mathrm{loc}+e^Y. ]

def fisk_rvs_numpy(
    c: float,
    size: int,
    loc: float = 0.0,
    scale: float = 1.0,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Sample Fisk(c, loc, scale) using NumPy only (inverse CDF)."""

    c = float(c)
    loc = float(loc)
    scale = float(scale)

    if c <= 0:
        raise ValueError("c must be > 0")
    if scale <= 0:
        raise ValueError("scale must be > 0")

    if rng is None:
        rng = np.random.default_rng()

    u = rng.random(size)

    # Avoid u=0 or u=1 exactly (would map to +/-inf)
    eps = np.finfo(float).eps
    u = np.clip(u, eps, 1.0 - eps)

    z = np.exp(_logit(u) / c)
    return loc + scale * z


# Monte Carlo validation (choose c>2 so mean/variance exist)
n = 200_000
samples_numpy = fisk_rvs_numpy(c0, size=n, loc=loc0, scale=scale0, rng=rng)

mc_mean = samples_numpy.mean()
mc_var = samples_numpy.var(ddof=0)

print("theory mean", m0["mean"], "MC", mc_mean)
print("theory var ", m0["var"], "MC", mc_var)

# Compare NumPy-only sampler to SciPy sampler (quick 2-sample KS test)
samples_scipy = stats.fisk.rvs(c0, loc=loc0, scale=scale0, size=n, random_state=rng)
ks = stats.ks_2samp(samples_numpy[:30_000], samples_scipy[:30_000])
ks

theory mean 2.13791866423119 MC 2.137557641110836
theory var  0.7145293838425228 MC 0.7035045300612408

KstestResult(statistic=0.009700000000000042, pvalue=0.11809773491898767, statistic_location=2.1769456563992517, statistic_sign=-1)

8) Visualization

We’ll visualize:

• the PDF with a Monte Carlo histogram overlay

• the CDF vs an empirical CDF

• the log-scale view ((\log X) should look logistic)

Because Fisk can be heavy-tailed, a log-x axis is often the clearest.

# PDF + histogram + CDF/ECDF
x_grid = np.logspace(-3, 2, 600) * scale0

# Empirical CDF (subsample for plotting)
sub = np.sort(samples_numpy[:10_000])
ecdf_y = np.arange(1, sub.size + 1) / sub.size

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=["PDF with Monte Carlo histogram", "CDF vs empirical CDF"],
)

# Histogram (density)
fig.add_trace(
    go.Histogram(
        x=samples_numpy,
        histnorm="probability density",
        nbinsx=120,
        name="samples",
        opacity=0.5,
    ),
    row=1,
    col=1,
)

# Theoretical PDF
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fisk_pdf(x_grid, c=c0, loc=loc0, scale=scale0),
        mode="lines",
        name="theory pdf",
        line=dict(width=2),
    ),
    row=1,
    col=1,
)

# Theoretical CDF
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=fisk_cdf(x_grid, c=c0, loc=loc0, scale=scale0),
        mode="lines",
        name="theory cdf",
        line=dict(width=2),
        showlegend=False,
    ),
    row=1,
    col=2,
)

# Empirical CDF
fig.add_trace(
    go.Scatter(
        x=sub,
        y=ecdf_y,
        mode="lines",
        name="empirical cdf",
        showlegend=False,
    ),
    row=1,
    col=2,
)

fig.update_xaxes(type="log", title_text="x", row=1, col=1)
fig.update_xaxes(type="log", title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)

fig.update_layout(height=420, title=f"Fisk(c={c0}, scale={scale0})")
fig.show()


# Log-scale view: log(X) should look logistic
y = np.log(samples_numpy - loc0)
mu_hat = np.median(y)
# logistic scale estimate via IQR: IQR = 2*s*log(3)
q25, q75 = np.quantile(y, [0.25, 0.75])
s_hat = (q75 - q25) / (2 * np.log(3))

# Compare to SciPy logistic on log-scale
logistic_dist = stats.logistic(loc=mu_hat, scale=s_hat)

y_grid = np.linspace(np.quantile(y, 0.01), np.quantile(y, 0.99), 500)

fig2 = go.Figure()
fig2.add_trace(
    go.Histogram(x=y, histnorm="probability density", nbinsx=120, name="log samples", opacity=0.6)
)
fig2.add_trace(go.Scatter(x=y_grid, y=logistic_dist.pdf(y_grid), mode="lines", name="logistic fit"))
fig2.update_layout(
    title="On the log-scale, Fisk becomes logistic",
    xaxis_title="y = log(x)",
    yaxis_title="density",
    height=380,
)
fig2.show()

9) SciPy Integration (scipy.stats.fisk)

SciPy parameterization:

stats.fisk(c, loc=0, scale=1)

• c is the shape.

• loc shifts the support to (loc, ∞).

• scale rescales (and in the 2-parameter form with loc=0, it equals the median).

Key methods:

• pdf(x), logpdf(x), cdf(x), ppf(q)

• rvs(size, random_state=...)

• fit(data, ...) for MLE

dist = stats.fisk(c0, loc=loc0, scale=scale0)

x_test = np.array([0.5, 1.0, 2.0, 4.0, 8.0])

pdf_scipy = dist.pdf(x_test)
cdf_scipy = dist.cdf(x_test)

pdf_ours = fisk_pdf(x_test, c=c0, loc=loc0, scale=scale0)
cdf_ours = fisk_cdf(x_test, c=c0, loc=loc0, scale=scale0)

print("pdf close?", np.allclose(pdf_scipy, pdf_ours))
print("cdf close?", np.allclose(cdf_scipy, cdf_ours))

# Fitting (MLE) with SciPy; fix loc if you know the support starts at 0.
c_fit, loc_fit, scale_fit = stats.fisk.fit(samples_numpy[:30_000], floc=0)
print("fit (c, loc, scale) =", (c_fit, loc_fit, scale_fit))

pdf close? True
cdf close? True

fit (c, loc, scale) = (4.998041440572122, 0, 2.007897677685366)

10) Statistical Use Cases

10.1 Hypothesis testing / goodness-of-fit

Typical workflow:

1. Fit the distribution parameters (MLE).

2. Check fit with diagnostics:

   • KS test (quick but less sensitive in the tails)

   • QQ/PP plots (especially important for tail fit)

   • tail-specific checks (log-log plots, exceedance probabilities)

You can also compare candidates (e.g., Fisk vs lognormal vs Weibull) using AIC.

10.2 Bayesian modeling

There is no conjugate prior for ((c,s)), but Bayesian inference is straightforward with MCMC/VI. A useful modeling trick is to work with (y_i = \log x_i):

[ y_i \sim \mathrm{Logistic}(\mu=\log s,,\sigma=1/c), ]

and add the Jacobian term if you write the model on (x) directly.

Below we’ll do a simple grid posterior approximation with priors on (\log c) and (\log s).

10.3 Generative modeling

Because sampling is cheap, Fisk is convenient for generating positive heavy-tailed synthetic data. For example, in an accelerated-failure-time style model, you can let the scale depend on covariates:

[ \log X = \beta^\top x + \varepsilon,\quad \varepsilon \sim \mathrm{Logistic}(0, 1/c). ]

# Hypothesis testing + model comparison (AIC) on synthetic data
true_c, true_scale = 2.5, 1.7
x_data = fisk_rvs_numpy(true_c, size=1500, scale=true_scale, rng=rng)

# Fit Fisk
c_hat, loc_hat, scale_hat = stats.fisk.fit(x_data, floc=0)
ll_fisk = float(np.sum(stats.fisk.logpdf(x_data, c_hat, loc=loc_hat, scale=scale_hat)))
aic_fisk = 2 * 2 - 2 * ll_fisk  # 2 params (c, scale) since loc fixed

# Fit lognormal
s_logn, loc_logn, scale_logn = stats.lognorm.fit(x_data, floc=0)
ll_logn = float(np.sum(stats.lognorm.logpdf(x_data, s_logn, loc=loc_logn, scale=scale_logn)))
aic_logn = 2 * 2 - 2 * ll_logn

# Fit Weibull
c_w, loc_w, scale_w = stats.weibull_min.fit(x_data, floc=0)
ll_w = float(np.sum(stats.weibull_min.logpdf(x_data, c_w, loc=loc_w, scale=scale_w)))
aic_w = 2 * 2 - 2 * ll_w

print("Fitted Fisk   (c, scale)", (c_hat, scale_hat), "AIC", aic_fisk)
print("Fitted lognorm(s, scale)", (s_logn, scale_logn), "AIC", aic_logn)
print("Fitted Weibull(c, scale)", (c_w, scale_w), "AIC", aic_w)

# KS test against fitted Fisk CDF (note: p-values are approximate after fitting)
dist_hat = stats.fisk(c_hat, loc=0, scale=scale_hat)
ks = stats.kstest(x_data, dist_hat.cdf)
ks

Fitted Fisk   (c, scale) (2.459344430482952, 1.6658920422291383) AIC 4812.347844437876
Fitted lognorm(s, scale) (0.7263480358234555, 1.661320104043914) AIC 4824.475216282033
Fitted Weibull(c, scale) (1.3345162113961129, 2.3881880476861284) AIC 5114.144353642654

KstestResult(statistic=0.01612416741026612, pvalue=0.8241228923230904, statistic_location=2.6508596780817344, statistic_sign=-1)

# Simple Bayesian grid posterior over (c, scale) with log-normal-ish priors
x_small = x_data[:200]
logx = np.log(x_small)

c_grid = np.linspace(0.8, 6.0, 120)
scale_grid = np.linspace(0.3, 4.0, 140)
log_scale_grid = np.log(scale_grid)

# Priors: log c ~ N(log 2, 0.7^2), log scale ~ N(log 2, 0.7^2)
mu0 = np.log(2.0)
s0 = 0.7

log_prior_c = stats.norm.logpdf(np.log(c_grid), loc=mu0, scale=s0)
log_prior_scale = stats.norm.logpdf(log_scale_grid, loc=mu0, scale=s0)

log_post = np.empty((c_grid.size, scale_grid.size), dtype=float)

for i, c in enumerate(c_grid):
    # Vectorized over scale_grid
    logz = logx[:, None] - log_scale_grid[None, :]
    a = c * logz

    # logpdf for each observation and scale
    logpdf = (
        np.log(c)
        - log_scale_grid[None, :]
        + (c - 1.0) * logz
        - 2.0 * _log1pexp(a)
    )

    loglik = logpdf.sum(axis=0)
    log_post[i, :] = loglik + log_prior_c[i] + log_prior_scale

# Normalize to probabilities
log_post -= np.max(log_post)
post = np.exp(log_post)
post /= post.sum()

# Posterior summaries
c_mean = float((post.sum(axis=1) * c_grid).sum())
scale_mean = float((post.sum(axis=0) * scale_grid).sum())

ij_map = np.unravel_index(np.argmax(post), post.shape)
c_map = float(c_grid[ij_map[0]])
scale_map = float(scale_grid[ij_map[1]])

print("posterior mean (c, scale) =", (c_mean, scale_mean))
print("MAP           (c, scale) =", (c_map, scale_map))

fig = go.Figure(
    data=go.Heatmap(
        x=scale_grid,
        y=c_grid,
        z=post,
        colorscale="Viridis",
        colorbar=dict(title="posterior"),
    )
)
fig.update_layout(
    title="Grid posterior p(c, scale | data) (loc fixed to 0)",
    xaxis_title="scale",
    yaxis_title="c",
    height=420,
)
fig.show()

posterior mean (c, scale) = (2.5669389607324895, 1.697934398969099)
MAP           (c, scale) = (2.547899159663866, 1.6841726618705037)

# Generative modeling sketch: scale depends on covariates (AFT-style)
n_gen = 1500
x_feat = rng.normal(size=n_gen)

beta0, beta1 = 0.2, 0.9
scale_x = np.exp(beta0 + beta1 * x_feat)  # positive, log-linear

u = rng.random(size=n_gen)
u = np.clip(u, np.finfo(float).eps, 1.0 - np.finfo(float).eps)

# Same shape c0, but varying scale per observation
x_gen = scale_x * np.exp(_logit(u) / c0)

fig = go.Figure()
fig.add_trace(
    go.Scatter(x=x_feat, y=x_gen, mode="markers", marker=dict(size=4, opacity=0.5))
)
fig.update_layout(
    title="Synthetic data: Fisk noise with covariate-dependent scale",
    xaxis_title="feature x",
    yaxis_title="generated positive outcome",
    height=380,
)
fig.show()

11) Pitfalls

• Invalid parameters: c <= 0 or scale <= 0 is not valid.

• Moment non-existence:

  • if c <= 1, the mean is infinite/undefined

  • if c <= 2, the variance is infinite

  • more generally, (\mathbb{E}[X^k]) exists only for k < c

• Parameterization confusion:

  • “Fisk” and “log-logistic” are the same distribution.

  • Many texts use (\alpha) (shape) and (\beta) (scale); SciPy uses c and scale.

  • On the log-scale, the logistic scale is 1/c.

• Numerical overflow:

  • direct computation of ((x/s)^c) can overflow for large x.

  • prefer logpdf, and compute the CDF via sigmoid(c * log(x/s)).

• Sampling edge cases:

  • inverse CDF uses logit(u); clip u away from 0 and 1.

• Fitting with loc free:

  • estimating loc changes the support boundary and can be unstable; fix loc when it is known (e.g., 0 for strictly positive data).

# Numerical stability demo: naive power vs log-space

def fisk_pdf_naive(x: np.ndarray, c: float, scale: float = 1.0) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    z = x / scale
    return (c / scale) * (z ** (c - 1.0)) / (1.0 + z**c) ** 2


x_big = np.array([1e2, 1e10, 1e50])
c_test = 5.0

with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
    naive = fisk_pdf_naive(x_big, c=c_test, scale=1.0)

stable = fisk_pdf(x_big, c=c_test, scale=1.0)
log_stable = fisk_logpdf(x_big, c=c_test, scale=1.0)

print("x", x_big)
print("naive  pdf", naive)
print("stable pdf", stable)
print("logpdf    ", log_stable)

x [1.e+02 1.e+10 1.e+50]
naive  pdf [0. 0. 0.]
stable pdf [0. 0. 0.]
logpdf     [ -26.0216 -136.5457 -689.1661]

12) Summary

• fisk is the log-logistic distribution: (\log X) is logistic.

• It models positive heavy-tailed data; the tail index is c.

• Moments exist only up to order < c (mean requires c>1, variance requires c>2).

• Sampling is easy via the inverse CDF: x = loc + scale * exp(logit(u)/c).

• SciPy provides a robust implementation: scipy.stats.fisk (pdf/cdf/rvs/fit).

References

• SciPy: scipy.stats.fisk

• Burr Type XII distribution (Fisk as a special case)

• Logistic distribution entropy and quantiles (for the log-scale view)